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Introduction. — Finding order in the apparent disorder 
of ocean motion is still an open problem, as even flows 
that in their Eulerian description are smooth present 
messy float trajectories [2]. Typical ocean structures are 
eddies and currents. Major currents, such as the Gulf 
stream or the Kuroshio, impact the Earth's climate be- 
cause of the heat they transfer. Eddies or rings are cer- 
tain robust, long-lived structures that may travel hun- 
dreds to thousand of kilometres, and persist for periods 
lasting from months to years. Understanding transport 
across these large scale structures is a challenging task, 
but one increasingly amenable to treatment since now 
data are becoming available [2 . Lagrangian tools provide 
a characterization of fluid flows. Underlying their de- 
scription is Poincare's idea of seeking geometrical struc- 
tures on the ocean surface (the phase portrait) that can 
be used to organize particles schematically by regions cor- 
responding to qualitatively different types of trajectories. 
For stationary flows the fixed points are key for describing 
the solutions geometrically. Fixed points may be classi- 
fied as hyperbolic or non-hyperbolic depending on their 
stability properties. Hyperbolic fixed points are responsi- 
ble for particle dispersion and non- hyperbolic fixed points 
are related to particle confinement. The interplay be- 
tween dispersion and confinement is an essential element 
of fluid transport processes. Stable and unstable mani- 
folds of hyperbolic fixed points divide the phase portraits 
in regions with qualitatively different types of trajectories 
since they are barriers to transport. This letter describes 
a new Lagrangian descriptor that for flows with a gen- 
eral time dependence realizes Poincare's idea of dividing 
a phase portrait in different regions that correspond to 
trajectories with qualitatively different behaviours. Our 
new instrument is based on a function which has been 
introduced in as a building block of a new definition 
of a Distinguished Trajectory (DT), which is a generaliza- 
tion of the concept of fixed point for aperiodically time 
dependent flows. Our function reflects, at the level of 
the phase portrait, relevant dynamical features of arbi- 
trary time dependent dynamical systems. Some of these 
features have not previously been detectable, thus, when 
applied to altimetric ocean data sets, M reveals the hid- 
den geometry of the ocean flow. The technique locates 
simultaneously hyperbolic and non-hyperbolic flow re- 
gions. Since it reveals singular features along the stable 
and unstable manifolds of the Distinguished Hyperbolic 
Trajectories (DHTs), it is useful as well to detect these 



invariant curves, that near the DHTs coincide with its 
stable and unstable directions. 

The function M. — The function we propose as a global 
Lagrangian descriptor (see [3 ), considers the system: 

x = v(x,t), x e M n , t e R (l) 

where v(x, t) is C r (r > 1) in x and continuous in t. 
Let x(t) denote a trajectory and denote its components 
in W 1 by (xi, x<i, x n ). For all initial conditions x* in 
an open set 8 G M n , at a given time t*, we define the 
function M(x*,f ) v , r : (B,i) ->> R for the system ([!]) as 
follows: 

M(x*,t*) v , r = y ( J2(dxi(t)/dt) 2 J dt (2) 

In our observational oceanographic flow, particle advec- 
tion occurs mainly in 2D (see [4 j), so n = 2 in 0. M is 
then the function that measures the Euclidean arc-length 
of the curve traced by a trajectory passing through x* at 
time t* on the plane (xi, X2). The trajectory is integrated 
from t* — r to t* + r. The function M depends on r and 
also on the vector field v. It is applicable to both time 
dependent and stationary flows. In the latter case it pro- 
vides a time independent partition of the phase portrait. 
For instance, for the unforced, undamped Duffing equa- 
tion a contour plot of M depicts the familiar stable and 
unstable manifolds of the fixed point located at the ori- 
gin. For time dependent flows the phase space partition 
provided by M is time dependent. 

The data set and the dynamical system. — The velocity 
data used in this work are geostrophic surface currents 
computed at CLS Int Corp (www.cls.fr) in the framework 
of the SURCOUF project |5]. The data span the whole 
Earth, at a resolution of 1/3 degrees at Equator, but 
we focus our results over a region through which passes 
the Kuroshio current, in selected days of May and June 
2003. Details on the data may be found in [6]. It is in- 
terpolated following methods described in [7J |H] , that use 
bicubic interpolation in space and Lagrange polynomials 
in time. Our coordinate system ((/>, /i), is related to the 
longitude and latitude (</>, A) by means of a transforma- 
tion \i = /i(A) (see details in [6[|7]). These variables are 
convenient for they distribute the data on a uniform grid. 
The equations of motion for ((/>, fi) are, 

# u(<j>,ii,t) dp v{<t>,^t) 

dt #cos(A(/i))' dt Ecos(A(/i)) 1 } 




FIG. 1: Evaluation of the function M over the Kuroshio current between longitudes 148°E-168°E and latitudes 30°N-41.5°N; 
a) and b) on May 2, 2003; c) and d) on June 3, 2003. Panels a) and c) take r = 15; Panels b) and d) take r = 30. 
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FIG. 2: The function M on May 2, 2003. a) r = 2; b) r = 15; c) the same as b) with a piece of stable manifold (black) and a 
piece of unstable manifold (green) of the DHT overlapping. 



where u and v represent the eastward and northward 
components of the altimetry surface velocity field re- 
spectively, and R is the radius of the Earth. The factor, 
l/cos(A(/i)), in the \i equation is an artefact of the coordi- 
nate transformation. The function M in Eq. Q is com- 
puted over the dynamical system (J3|. Thus the length 
of the trajectory is measured on the (0, \i) plane. The 
system expressed in Eq. ([3| is not exact, as it is subject 
to errors coming from the measured velocity fields, the 
sort of interpolation used, etc. However Eq. ([3| is used 
for evaluating M, a function that contains Lagrangian 
information. In the literature [9], has been studied the 
robustness of the Lagrangian structures under errors in- 
duced in the vector field satisfying certain conditions. We 
have assessed the reliability of M by computing it with 
several interpolation schemes. 

Results. — We demonstrate that the function denned in 
Eq. ([2| gives a global dynamic picture of oceanic flows 
since it detects simultaneously invariant manifolds, hy- 
perbolic, and no n- hyperbolic flow regions. It synthesizes 



information more efficiently than for instance spaghetti 
diagrams. These represent paths over time of messy tra- 
jectories but they do not communicate information about 
regions in which particle evolutions are qualitatively dif- 
ferent, and one cannot get much intuition from them. 
Figure [I] displays the function M for medium and large 
r on selected days of May and June 2003 along the me- 
andering Kuroshio current. Maximum values of M are 
in red, while dark blue indicates minima. The depen- 
dence of M on time is obvious for this highly aperiodic 
flow, since representations for different days have differ- 
ent structures. In the figure the organizing centres are 
visible at a glance. These key points are the minima 
of M, and as discussed in [3] they are related either to 
hyperbolic or non-hyperbolic distinguished trajectories. 
Singular features of M forming lines are easily discerned, 
both in Figs. [I] and [5] Fig. |2Jd) shows their intersection 
at a hyperbolic minimum at longitude ~ 157.1° and lati- 
tude ~ 35.63°. The time evolution of this point has been 
characterized in [6] as a DHT. Singular lines are identi- 
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FIG. 3: The function M evaluated over the inner core of an eddy on May 2, 2003. a) r — 15; b) r = 30; c) r = 72. 



fied as manifolds since they are advected by the flow and 
are asymptotically obtained from small segments aligned 
with the stable and unstable subspaces of the DHT. Fig. 

shows the overlapping of M with the stable and un- 
stable manifolds computed with the technique used in 
[4]. This confirms the coincidence of the lines with the 
manifolds. Why should stable and unstable manifolds 
be traced out by singular features of M? M measures 
the lengths of curves traced by trajectories on the phase 
space, so it is expected it will change abruptly at the 
boundaries of regions comprising trajectories with quali- 
tatively different evolutions, since this is exactly what the 
stable and unstable manifolds separate. Convergence of 
the structure of M towards these singular lines requires 
a large enough r value. For instance for r = 2 the ap- 
pearance of M in Fig. |2^i) is rather simple, almost with- 
out structure and resembling that of Eulerian currents, 
while sharp lines in Fig. ^p) and c) require the use of 
r = 15. The structure of M becomes more and more re- 
fined for larger r values as confirmed by panels [TJd) and 
d), obtained for r = 30. This is justified because M re- 
flects the history of initial conditions on open sets, and 
in highly chaotic systems this history is expected to be 
more complex for longer time intervals. The evaluation 
of M in large oceanic areas as shown in Fig. [T] reveals 
recognizable phase portraits similar to those of the cat's 
eyes of the forced pendulum (in panel a) upper left), or 
the forced Duffing equation (see panel c) at the lower 
right). The ocean surface resembles a patchwork of inter- 
connected dynamical systems from which the complexity 
of possible particle routes is visible. 

In figures [l] and |2j apart from the minima of the func- 
tion M at the intersections of singular lines, related to the 
hyperbolic DT, there are apparent minima at the eddy 
centres. In the work by Madrid and Mancho [3 these 
have been related to non- hyperbolic DT (DET), which 
are eddy-like structures, of great interest to oceanogra- 
phers. The Lagrangian description of eddies, such as that 
shown in Fig. [2] reveals the existence of an inner core, 
which is robust and rather impermeable to stirring and 
an outer ring, where the interchange with the media is 
understood in terms of lobe dynamics (see [E]). We ana- 
lyze how the function M reflects to what extent the inner 
core of Fig. [2] is impermeable to mixing. Figure [3] dis- 
plays contour plots of M on t =May 2, 2003 for several r. 
In Fig. pa) it is observed that for r = 15 days the interior 



of the eddy has a minimum which is locally smooth. This 
means that in the range (t — r, t + r) trajectories in this 
neighbourhood outline similar paths, and for this reason 
the function M does not change sharply (i.e does not 
have singular features). Smoothness of M implies that 
for these initial conditions it does not perceive nearby 
hyperbolic regions for (t — r, t + r). Hyperbolic trajecto- 
ries are the ones responsible for dispersion and it is just 
these trajectories that may induce sharp changes in M. 
In Figures ^jp) and c) for larger r values (i. e, r = 30 
and 72 respectively) the interior of the eddy becomes less 
and less smooth, for in the range (t — r, t + r) trajectories 
placed at the interior core either were dispersed in the 
past or will disperse in the future. In fact, in Figure [3 3), 
the interior of the core is completely foliated by singular 
features associated either to stable or unstable manifolds 
of nearby hyperbolic trajectories. So, the value at which 
M starts losing smoothness, e.g. 2r = 60, is a good indi- 
cator of the maximum time for confinement of particles 
in the inner core. The minimum of M on the elliptic re- 
gion does not converge with r, and this is the condition 
required for finding DT. Similarly to what is described in 
[3], DET have not been found in highly aperiodic flows. 
Figure [3^) displays in black line with a computed un- 
stable manifold which overlaps on the contour plot of 
M. Again there is observed a coincidence of the singular 
features of M with the manifold. However the foliated 
structure of M is much richer than that provided by the 
manifold. The reason is that the manifold has been com- 
puted from the one DHT recognizable in Fig. 2b) and 
c), while M displays all stable and unstable manifolds 
from all possible DHTs in the neighborhood of the eddy, 
without need for identifying DHTs a priori, as required 
by the manifold algorithm (see [4 ). Thus M provides a 
complete partition of the phase portrait, while the direct 
computation of a manifold of a DHT does not. 

The Lagrangian method using M has several advan- 
tages over other methods based on finite time versions 
of Lyapunov exponents (LE) such as FTLE or FSLE 
[TTj H2] . LE techniques provide information on the lin- 
earised flow along trajectories and their focus is on hyper- 
bolic regions. Ridges of FTLE and FSLE fields represent 
manifolds as reported for instance in [TTJ [12] . Figure [4] 
confirms this point. In it there is displayed the same 
eddy of Fig. |3J but for r = 50. The sharpest ridges 
of the FTLE field represented in Fig. Hp) and c) are 
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associated to the stable (forward FTLE) and unstable 
(backward FTLE) manifolds respectively, in close corre- 
spondence to the singular features of M displayed in Fig. 
^l). An obvious difference between these two representa- 
tions is that M contains the information on both stable 
and unstable manifolds in the same picture, while FTLE 
splits it in two diagrams. Regarding other features pro- 
vided by FTLE and M there is not strict agreement. In 
the regions centered at longitude ~ 156.1° and latitude 
~ 36.8° and at longitude ~ 156.3° and latitude ~ 36.5°, 
the structure of M is smooth and eddy-like, which as 
explained above, indicates that particles in that area do 
not disperse. The same areas in Fig. [4] b) display a 
striped pattern suggesting that particles disperse in this 
time interval. Numerical integrations of trajectories in 
these regions confirm that they stay close to each other 
in the interval (t — r, t + t). Non-sharp striped structures 
in Fig. [4] b) raise then the questions of the kind of in- 
formation they provide and if there is an upper bound 
on r for this vector field, beyond which the validity of 
some of the structures provided by the FTLE cannot 
be confirmed. In [13 it is reported that for particular 
transient flows FTLE may develop 'ghost' structures, al- 
though a detailed discussion on this is beyond the scope 
of our letter. A further difference is that FTLE require 
certain assumptions on the vector field (see [H]), while 
by construction the function M is defined for a general 
time dependent vector field and there is required only 
existence and uniqueness of solutions for the system ([3]). 



Computationally the evaluation of the function M is 
cheaper than that of LE. In geophysical flows, both M 
and LE require the performance of a large number of in- 
tegrations on a dynamical system such as (|3|, where the 
vector fields are interpolated over a finite space- time grid. 
Interpolations make computations expensive and saving 
these at each time step is a convenient feature. Evalua- 
tion of the function M fits this criterion better than LE 
as each point in the phase space requires of just one inte- 
gration forwards and backwards in time. LE techniques 
however require more interpolations of the velocity field 
at each time step, either because they evaluate a sepa- 
ration rate among several trajectories, or because they 
compute the linearised flow around each trajectory. 

The Lagrangian descriptor M locates special organiz- 
ing trajectories called DT as reported in [3]. However 
although there are references suggesting the ability of LE 
to achieve this goal [14], there are no published studies 
where this is discussed in detail. 

In conclusion, this work demonstrates the efficiency of 
a new Lagrangian descriptor M, for identifying the essen- 
tial dynamical features of general time dependent flows. 
M is a promising tool for the estimation of transport in 
realistic flows, as recent articles [15] have confirmed. 
Pursuing further the discussion on how it contributes to 
transport diagnosis, is beyond the scope of this letter. 
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